import matplotlib.pyplot as plt
import numpy as np
#函数1
def func1(x,y):
    return (5*y-6*np.exp(-x))
#函数2
def func2(x,y):
    return (x+y)*np.sin(x*y)
#欧拉法
def ola(x,y0,h):
    y=np.zeros(len(x))
    #初值
    y[0]=y0
    for i in range(len(x)-1):
        #根据上一个值求下一个值
        y[i+1]=y[i]+h*func2(x[i],y[i])
        
    return y
#三阶龙格-库塔法
def lunge_3(x,y0,h):
    y=np.zeros(len(x))
    y[0]=y0
    k1,k2,k3=0,0,0
    for i in range(len(x)-1):
        k1=func2(x[i],y[i])
        k2=func2(x[i]+h/2,y[i]+h*k1/2)
        k3=func2(x[i]+h,y[i]-h*k1+2*h*k2)
        y[i+1]=y[i]+(k1+4*k2+k3)*h/6

    return y
#四阶龙格-库塔法
def lunge_4(x,y0,h):
    y=np.zeros(len(x))
    y[0]=y0
    k1,k2,k3,k4=0,0,0,0
    for i in range(len(x)-1):
        k1=func2(x[i],y[i])
        k2=func2(x[i]+h/2,y[i]+h*k1/2)
        k3=func2(x[i]+h/2,y[i]+h*k2/2)
        k4=func2(x[i]+h,y[i]+h*k3)
        y[i+1]=y[i]+(k1+2*k2+2*k3+k4)*h/6
        
    return y
#问题1
fx0=1.0+1.0e-13
x1=np.arange(0,6,0.01)
x2=np.arange(0,6,0.001)
x3=np.arange(0,6,0.0001)
#问题2
fx0_=5
x4=np.arange(0,6,0.2)
x5=np.arange(0,6,0.1)
x6=np.arange(0,6,0.01)

#y=ola(x6,fx0_,0.01)
#y=lunge_3(x6,fx0_,0.01)
y=lunge_4(x6,fx0_,0.01)
print("求出的数值解为：\n",y)
plt.plot(x6,y)
plt.xlabel("x")
plt.ylabel("y")
#plt.title("(b)欧拉法,h=0.01")
#plt.title("(b)三阶龙格-库塔法,h=0.01")
plt.title("(b)四阶龙格-库塔法,h=0.01")
#使标题输出为中文
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.show()